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^ \ We present an algorithm to analyze numerically the bounce solution of first-order 
, ^3 ■ 

phase transitions. Our approach is well suited to treat phase transitions with several 

(N 

CN . fields. The algorithm consists of two parts. In the first part the bounce solution 

\ without damping is determined, in which case energy is conserved. In the second 



part the continuation to the physically relevant case with damping is performed. 
The presented approach is numerically stable and easily implemented. 
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Ph! I. INTRODUCTION 

Oh: 

• The problem of calculating properties of a first-order phase transition from a false vacuum 

- ^ , to an energetically lower vacuum appears in different contexts of physics, as for example con- 

' densed matter physics , particle physics 0, 0] , and cosmology P, 0, Q] . In a cosmo logical 



context, one of the interesting features of a first-order phase transition is that it can proceed 
by bubble nucleation. The subsequently expanding bubbles can be utilized as a source that 
drives the hot plasma in the early Universe out of equilibrium. This is in turn one of the 
Sakharov conditions, the requirements for a viable baryogenesis mechanism. Employing the 
electroweak phase transition, the corresponding baryogenesis mechanism was suggested in 

genesis. Even though this scenario of the generation 
of the baryon asymmetry of the Universe (BAU) is excluded in the Standard Model (SM) 
because of a lack of CP violation and the cross-over nature of the phase transition, this 
possibility is still open in the Minimal Supersymmetric Standard Model (MSSM) jol. Ilfll. Illl| . 
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An essential input parameter of the determination of the produced BAU is the wall shape 
of the expanding bubbles. While in the SM the phase transition is described by a single vev 
of the Higgs field v, there are already two dynamical fields during the phase transition in 
the MSSM, namely two vevs, vi and V2, or equivalently, v = a/w i + f | and tan(/3) = V2/V1. 
Notice that the initial value of (3 cannot be inferred without studying the dynamics of the 
phase transition, since v vanishes in the symmetric phase. This is especially important, 
because some CP violating sources are proportional to the change of (3 during the phase 
transition. 

Examining extensions of the MSSM, one has eventually to deal with more than two scalar 
fields. TK. . the case tKe MSSM extended a„ add,t,o.a. si„g.et UUU, 

which provides the interesting feature of transitional CP violation Further examples 
are provided by models that equip this singlet with a f/(l)' gauge symmetry or by the 



SM with two Higgs doublets 
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ins rates in particle physics 
|. However to calculate 



The theoretical basis for the determination of the tunne 
and the explicit solution for one scalar field was given in 
the properties in the case of several scalar fields is analytically not accessible, and even 
numerically a non-trivial task. We present a numerical solution to this problem. 

II. THE PROBLEM 

Before the numerical treatment of the multi field case is discussed, we will briefly review 
the results of 0] and clarify our notation. Given the Lagrangian 

c = ^d,<j,d^<j>-v {<!>), (1) 

the transition probability per unit time and unit volume from a false vacuum to an energet- 
ically lower vacuum is, in the semi-classical approximation, given by 
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x(l + 0(n)). (2) 



det[-d^ + V"( 

The primed determinant indicates that the zero-modes, corresponding to translation in- 
variance, have been omitted, leading to the factor S[(f)]'^ in the formula. The function is 
the classical Euclidean solution of the equations of motion derived from the Lagrangian in 
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Eq. (PJ) with inverted potential and the boundary conditions 

hm 0(r, x) = 0_, (3) 

where 0_ denotes the false vacuum and 

dr(j)iO, x) = 0. (4) 

In addition the solution should be non-trivial in the sense that 0(0, x) = (ph is close to 
the new vacuum. The corresponding solution is usually called bounce configuration and a 
typical potential is illustrated in Fig. [H 
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FIG. 1: The potential for a theory with a false vacuum. 

The present work aims at determining numerically the bounce configuration in theories 
with several scalar fields. Thus, the parameters that are inferred about the phase transition 
rely on the semi-classical ex pan sion in Eq. 0, which can break down, e.g. for very weak 
first-order phase transitions [l9|. 

In addition, the potential is treated as static and it is assumed that the phase transition 
occurs at some fixed temperature. Issues concerning the timing of the phase transition can 
be decided on by comparing the nucleation rate with the corresponding cooling time scales, 
such as the expansion rate of the Universe. This includes, for instance, the questions to 
know if the vacuum resides in a metastable state until today or if the phase transition is 
performed before the Universe cools down to a temperature where only a second-order phase 
transition is possible 
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A crucial question in this context, relevant to electroweak baryogenesis, is to see if the 
phase transition is strong enough to avoid the washout of the generated baryon asymme- 
try [2^. Fortunately, it turns out that the perturbative effective potential even slightly 
underestimates the strength of the phase transition in the MSSM 0, 

Finally, for certain physical problems the knowledge of the prefactor in the decay rate of 
Eq. Q might be necessary, which we do not attempt to achieve. For a recent development 
and a review on the topic, see Refs. Q, Q- 

Since all these considerations are model- 
dependent, we will focus in the following on the numerical determination of the bounce 
solution. 

The solution for one scalar field in the thin-wall regime can be analytically determined 
by using the saddle-point method 0, 14 1 . Using the fact that there exists a solution that is 
invariant under Euclidean four- dimensional rotations, and hence that only depends on the 
combination 



P 



the corresponding action and equation of motion (EoM) read 
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^[0(p)] = 27rM dpp- 







dp 



and 



dp^ p dp dcj) 



(5) 



(6) 



(7) 



Numerically, in the case with one field, this solution can be easily found by the 
"overshooting-/undershooting" method. Varying the initial point 0{, and testing if the 
configuration after integration of the EoM over- or undershoots the position 0_, the bounce 
point 0b and the corresponding configuration can be determined. Clearly this method can- 
not be generalized to the multi field case, since starting at 0+ there is no viable strategy to 
find the bounce point 0;,. 

For tunneling in a thermal system, the 0(4)symmetry is replaced by an 0(3) symmetry 



and the required solution is constant in time 
read 



symme 

[la. 



The corresponding action and EoM 



5[0(p)]=47r/ dpp 



dp 



+ V{<P,T) 



and 



f± + 2rf0 ^ —Vict) T) 
dp"^ p dp d(j) ' ' 



(8) 



(9) 
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and the "overshooting-/undershooting" method can be utihzed in a similar fashion. Gen- 
erahzing these two cases and allowing for n scalar fields, the corresponding action yields 

2 



~sMp)] = / ^pp^"-'^ 



1 

2 \dp 



+ ^(0) 



(10) 



where a = 3 (a = 4) in the temperature (vacuum) case, denotes a vector with n entries, 
and a constant prefactor of the action is neglected. The n equations of motion read 

^''^ + = VK(0) (11) 



dp^ p dp 

and can be interpreted as a classical particle moving in the tt,- dimensional inverted potential 
with a time-dependent damping term. In the thin-wall regime, where the maxima are almost 
degenerate, 1^(0+) ~ y(0_), one can determine the action analogously to that in Ref. 
as 

7ft AT/ 

^a[p] = -^^ + p("-')^i, (12) 
a 

where 5V denotes the potential difference between the two maxima and p the radius of the 
bubble. Minimization of the action leads to the radius of the bubble 

P=(«-l)^, (13) 
and yields the following estimate of the action in the thin-wall regime 



SMp)] = ^ Jy.-. (14) 



The problem is now to determine the bounce solution apart from the thin-wall regime. 
Even with only two fields, it is a challenging task to find numerically the position (pb where 
the fields bounce. One possibility is, starting at 0_, to "scan" into different directions. For 
two fields this seems to be feasible, but it clearly fails for more than two fields. 

Alternatively, one could try to obtain the bounce solution by minimizing the action. 
However, it was already pointed out by Coleman that the bounce is not a minimum of the 
action but a saddle point [26]. In the following paragraphs we will briefly review former 
techniques to solve this problem, before we discuss our approach in Sections IIIII to El 



The first improvement was made in Ref. |22|, where a tunneling amplitude in supergravity 
is discussed. The authors identified the direction of the negative eigenvalue of the action 
to be the scaling mode of the field (p. Accordingly, they minimized the action in several 
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iterations and maximized it with respect to scaling, simultaneously. However, to determine 
the directions orthogonal to scaling is numerically a difficult task. The authors reported 
errors of 3-5% in the solution, which is due to the fact that the system converges to a point 
where scaling and minimizing balance, and not to the real bounce solution. 

This weakness was criticized in Ref. j^l, where another paradigm was followed, the so- 
called improved-action approach. The direction of the negative eigenvalue was lifted in the 
vicinity of the bounce by using a generalized version of the viral theorem. To implement 
this procedure successfully, one needs a configuration that is already close to the bounce. 

Further approaches are based on minimizing an action that consists of the squared equa- 
tions of motion This method also depends on a good anticipation of the bounce. In 
addition, it was reported in Ref. Q| that in certain cases the implementation is numerically 
not stable in the sense that convergence to the bounce is not guaranteed and an oscillatory 
behavior and unphysical solutions can evolve. 

The last two methods are not well suited to problems with more than two fields, since a 
good anticipation of the bounce is in this case usually not possible. 

Another method based on performing a certain order of gradient descent and ascent steps 
was applied in Ref. |31i]. In this algorithm, the ascent steps ensure that the configuration 
approaches a local extremum and not just a local minimum. Additionally, a second algorithm 
based on a combination of the shooting and the gradient-descent method was presented in 
this work. This second algorithm is slightly more technically involved than our proposal, 
since it uses spline interpolation. These two algorithms have a similar range of validity 
as our approach, but they are slower, so it is more numerically expensive to achieve high 
accuracy. 

To find the bounce solution, we split the problem into two parts. First, we determine the 
solution without damping (a = 1). Subsequently, we perform a smooth transition to the 
damped case of interest (a = 3 or a = 4). Neglecting the "damping" term is essential in 
the first part, since the corresponding approach is based on the assumption that energy is 
conserved and the bounce point 0b is on the submanifold with potential equal to V{(f)-) (as 
already indicated in Fig. H}. 
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III. IMPROVED POTENTIAL APPROACH 

In this section we will present our approach to find the bounce solution in the undamped 
case, a = 1. 

First, notice that, even in the undamped case, the initial problem remains and the bounce 
configuration is only a saddle point of the action. Intuitively, this can be understood as 
follows. The undamped case can be interpreted as a particle moving freely in a potential. 
To minimize the action, the system tries to avoid valleys, and hence moves fast in them. If 
the system has a long time to evolve, the global minimum of the action is given by a field 
configuration that shoots up to the maximum of the inverted potential 0+ and stays there 
as long as possible, only restricted by the boundary conditions. Thus, the bounce solution 
cannot be a global minimum of the action. 

Furthermore, close to the bounce solution that fulfills the boundary conditions: 

0(0) = 0,, 0(oo) = 0_, (15) 

the action can be reduced by increasing 9^0(0) in the direction of maximal slope (towards 
0_i_). Hence the bounce configuration constitutes a saddle point of the action, even though 
the "damping" term is neglected. Notice that in the case of degenerate minima the negative 
eigenvalue corresponding to the instability of the saddle point [22| becomes zero and refiects 
the invariance under translation in the coordinate p. The solution of the degenerate case 
without damping can be easily obtained by minimization of the action, and the subtlety of 
instabilities appears only if the minima are non-degenerate. 

Our starting point for the determination of the bounce solution is the configuration 
restricted by the following boundary conditions at the maxima of the inverted potential: 

0(0) = 0+, 0(T) = 0_, (16) 

where T is a time scale that is sent to oo iteratively. Given that 0+ is the global maximum 
of the inverted potential, there exists a solution with these boundary conditions that indeed 
is a minimum of the action in Eq. (fTIHl and can therefore be found by simply discretizing 
the field configuration and minimizing the action using the Newton method. 

In the case of non-degenerate minima the resulting configuration is not the bounce. Our 
goal is to modify the potential in such a way that the minimization of the standard action 
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with modified potential leads first to a good prediction of (pi, and finally to the bounce 
configuration. With respect to the original problem, only changes in the potential are 
performed, and an interpretation in terms of a classically moving particle is viable throughout 
the present discussion. 

In a first step the potential is changed by performing the following transformation: 



t,.W^^W-^(^-)+J'^-W (17) 
2 V 4 

For finite e, this modifies the potential in a differentiable way and in the limit e — » 

transforms the potential in such a way that it is unchanged (besides a constant) for V ^ 
V{(j)-) and develops a constant plateau in the region with V ^ V{(j)J). 
In the simultaneous limits 

e 0, T ^ oo, (18) 

the minimum of the modified action will determine the bounce solution. Here, the uniqueness 
of the minimum of the action is assumed, even though in some cases this is not valid j^^l- 
To show that this minimum determines the bounce configuration, the configuration that 
corresponds to the minimum of the modified action is split into two parts. The first part 
consists of the path in the plateau, while the second part is outside the plateau, and both 
are connected by a point 0e = 0(te) that is defined by V{(f)e) = V{4'-)- The path in the 
plateau does not contribute to the action in the limit te ^ oo, e ^ 0, and hence does not 
restrict the point (pe or the path outside. In the limit e ^ 0, the path outside the plateau 
fulfills the Euler-Lagrange equations in the original potential V{(j)). Using the fact that, 
in this limit, U^{(j)J) and f/e(0e) both approach the global maximum U^{(j)^) of the inverted 
potential, leads to the fact that 

dp<P{Q^O, 9p0(T)^O (19) 

in the limit (T — tg) oo. Hence, the path outside the plateau defines a bounce solution in 
this limit according to Eqs. © and (jH). 

For finite e and T, the field configuration will move very slowly in the plateau from 0+ 
to 0e (a point at a distance 0(e, 1/T) to (pb) and then move from 0e to 0_, much as the 
bounce solution. This way, a fairly good estimate of the bounce configuration and of (pb 
is obtained even for finite e and T. On the other hand, this configuration still misses one 
typical feature of the bounce, namely the exponential behavior of close to 0_. Of course. 
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this feature appears in the hmit of Eq. ()18|). but the numerical stabihty and performance 
of the algorithm can be improved by an additional change of the potential. The reason for 
the absence of the exponential behavior is that, since the global minimum of the potential 
is still close to 0+, the configuration spends most of the time in the plateau. To avoid this 
effect, another term can be added to the potential, of the form 

f/,(0) = f/,(0) + A[/,(0) (20) 

This additional contribution is designed to lift the potential so that 

< - < 0(e) (22) 

and the corresponding solution shows the exponential behavior close to (j)-. This lifting 
is in principle not necessary, but it improves the performance of the algorithm drastically. 
Initially, a reasonable choice for e is 



-init 



r(0+)-V(0_), (23) 



while for T a time that approximately corresponds to 20 times the thickness of the wall is 
appropriate. The wall thickness can be anticipated by determining the local minimum of 
the inverted potential and |0_ — <t>+\ 

^init = 20 X 1^-"^;! _ ■ (24) 

By the procedure described above one obtains a good estimate 0e of the bounce point 
0{, and a configuration close to the desired bounce configuration. The convergence can be 
further improved by truncating the part of the configuration in the plateau. The boundary 
conditions are accordingly changed to 

0(0) = 0e, 0(T) = 0_, (25) 

where 0e initially denotes the position that corresponds to the point with V{(j)e) = ^(0-) 
on the field configuration that was obtained by taking the previously described steps. In 
the next minimization procedure, (pe is allowed to move freely on the submanifold with 
y{4>e) = ^(0-) SO that 0e ^ 0fe IS eusured. While minimizing the action with the potential 
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FIG. 2: The bounce configuration in the one-dimensional degenerate case. The deviations from 
the exact analytic result are O(10~^). 

Ue{4'), e can be iteratively sent to zero until the numerical uncertainties inhibit further 
improvement of the bounce configuration. 

To truncate the configuration in the plateau is helpful for two reasons. First, the energy is 
conserved and is exactly V {</)-) for the bounce solution. The motion in the plateau from 0e 
to 0+ contributes at least (0e — 0+)^/2T^ to the kinetic energy, which requires very large T 
for accurate results. This would lead to large lattices, provided the physical lattice spacing 
is kept constant. This effect is avoided by starting at 0e- Secondly, since U^{(j)) is very fiat 
in the plateau for e ^ 0, the convergence of the Newton method is very poor in this region. 
Starting at 0e, the field never enters the plateau. 

Even though the changes outside the plateau are only marginal for e — 0, it is clear, from 
classical reasoning, that this procedure stabilizes the unstable mode of the original action. 

IV. NUMERICAL RESULTS OF THE UNDAMPED CASE 
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As a first example, the degenerate case is demonstrated. Our approach is applied to the 
following potential with one scalar field and degenerate minima 

\/(0) = 2(0-l)V. (26) 
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FIG. 3: Potential V and kinetic energy E-^^^ of the one-dimensional configuration. The total 
energy E = E-^^^ — V is conserved up to 0.1%. 

The resulting configuration and the deviation from the exact analytic result are plotted in 
Fig. 121 The deviation from the exact solution is O(10~^) and mostly due to the finite lattice 
spacing. Another test is provided by the kinetic energy that should vanish in the limits 
p — ^ and p —>■ T. This was not enforced by the boundary conditions since the start and 
end-point of the configuration are used as boundary conditions of the discretized action, as 
given in Eq. ()16|) or Eq. ()25|) . In addition one can analyze the conservation of energy. These 
tests are displayed in Fig. El 

In the second example, a non-degenerate potential for the multi scalar case is analyzed, 
namely 

= 16(01 - + 2 02 + ci 01 + C2 020l(01 - 1). (27) 

The third term dissolves the degeneracy of the minima and is chosen to be rather small 
(ci = —0.1). The last term gives a non-trivial dynamics in the 02 direction (c2 = 8.0). The 
numerical results are depicted in Figs. and O The results show the expected exponential 
behavior close to 0_ and a parabolic shape near 0b . The energy is again conserved up to 
the per mille level. Especially the parabolic shape with dp(j) = at the bounce point 0^ in 
Fig. Hlis a clear indication of the high precision of our numerical approach, since because of 
the lift of the potential in Eq. (jUJ), the kinetic energy is there 0(e). 

Note that the shape of the field 02 is not of the usual hyperbolic tangent form, and 
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FIG. 4: The bounce configuration in the two-dimensional non-degenerate case. The exponentially 
decreasing tail [0, 250] has been cut off. 
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FIG. 5: Potential V and kinetic energy E^^^ of the two-dimensional configuration. The total 
energy E = E^^^^ — V is conserved up to 0.1%. The exponentially decreasing tail [0, 250] has again 
been cut off. 

not even monotonically increasing; hence one cannot obtain a good estimate of the bounce 
configuration by minimizing the action with respect to a specific ansatz, as was done in 
Ref. |30|. 



*l>2 
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V. CONTINUATION TO THE DAMPED CASE 

To perform the continuation from the undamped case (a = 1) to the damped case, a 
numerical method to solve the discretized version of Eq. (jllj) with the boundary conditions 
as given in Eqs. (jHl) and (jlj) is needed. In the present work the following discretized form of 
the EoM has been used: 



bi+i - 2(j)i + a -I (t)i+i - 0i_i 



VF(0,), ^G[l,n-1], (28) 



cip2 {i + /\i)dp 2 dp 

with the boundary conditions 

00 = 01, 0n = 0-, (29) 

where 0, constitutes the one-dimensional lattice {i G [0,^]). This discretization has the 
virtue to be derived from the energy functional 

i/ = E'M« + E('^ + o^^^^^ (30) 

i i 

and reproduces the continuum energy with an error of 0{dp^). 

The offset parameter Ai has been introduced to improve the performance of the algorithm, 
since an initial value of Az = can easily lead to a pathological behavior close to z = 0. 
Finally, this parameter has to be sent to zero. 



One way to obtain a solution of t 



of motion, as was done in Refs. 
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lis system, is to use the method of squared equations 



30|. However, since the algorithm in the last sections 
provides fairly good estimates (pi of the solution, we pursue another possibility. Using the 
Taylor expansion of the gradient field 

Vy(0O ^ Vy(0.) + (0. - 0.) ■ V{VV{4>^)), (31) 

we can linearize the system of equations and iteratively solve the equations by inversion 
of the linear system. Since the matrix corresponding to this system of equations is band 
diagonal, its inversion is easily performed. The original system of equations in Eq. (PHj) is 
then solved iteratively by using the linearized system. Notice that the original potential is 
used and not the modified potential that has been utilized in the minimization of the action 
in the last section. 

This method has been applied to the non-degenerate example of the last section and 
the results are briefly displayed in the following. Starting point of the continuation are the 



14 




FIG. 6: The fields 0i and (j)2 during the continuation for several values of a and Ai = 300. 




FIG. 7: Contour plot of the potential and several trajectories. 

damping parameter a = 1, an offset Ai = 300, a time scale dt = 0.0158, and a lattice of 
size n = 400. As initial configuration, the one generated by the improved-potential method 
is used, as discussed in the last section. 

The parameter a is increased by steps of size 6a = 5 x 10"'^ and the linear system in 
Eq. ()28j) is solved iteratively for each parameter set and some obtained configurations are 
plotted in Fig.lHl The contour of the potential and several trajectories are depicted in Fig. [7| 



As soon as the size of the bubble 

p={a-l)^-Azdt, (32) 
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lattice sites 


^3 


ASs 


100 


725.819 




200 


726.957 


1.137 


400 


727.239 


0.282 


800 


727.309 


0.070 


1600 


727.326 


0.017 



TABLE I: The action 5*3 for different numbers of lattice sites. The change in the action AS3 
behaves as 0{6p'^). The error of the largest lattice can be estimated to be 5 x 10~^ what corresponds 
to a relative error of 7 x 10~^. 

is larger than the wall thickness, an increase in a or a decrease in Ai results almost only in 
a translation of the configuration. If this translation is anticipated, a faster continuation in 
the parameters a and Ai becomes possible. 

Since solving the linearized system of equations is not very time consuming, the accuracy 
of the result is in principle only limited by the finite size of the lattice and the accuracy of 
the numerical representation. This can be demonstrated by increasing the number of lattice 
sites and the results are summarized in Table HI The relative error in the action behaves 
as expected as 0{6p^) and reaches approximately 7 x 10"^ for the largest lattice that was 
used. Errors of 0{Sp^) can be obtained by extrapolation of the action to vanishing 6 p. In 
comparison, using the method of squared EoM or one of the algorithms of Ref. js^ is much 
more time consuming and leads to inferior accuracy. 



VI. CONCLUSION 



We presented a numerically stable and easily implemented algorithm to determine prop- 
erties of phase transitions with several scalar fields. We demonstrated our approach in the 
simplest one-dimensional and in a non-trivial two-dimensional example. The approach we 
suggest converges with high precision to the bounce configuration. During the minimiza- 
tion neither oscillatory nor unphysical solutions have been encountered, in contrast to the 
approach followed in Ref. js^. 

Our approach was based on a procedure consisting of two steps. First, we determined the 
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bounce solution in the undamped case, using a novel algorithm. Secondly, a smooth contin- 
uation was performed to the physical case, turning on the damping term in the EoM. Note 
that the sole purpose of the undamped solution is to serve as an initial configuration for the 
continuation. In contrast to the configuration without damping, the damped configuration 
usually shares more characteristics with the typical solution in the degenerate case, namely 
a bubble size much larger than the bubble wall thickness and an exponential behavior close 
to both maxima. Hence the action determined in the thin-wall regime as given in Eq. ()14|) 
is typically more accurate than using the configuration of the undamped case in the action. 

The main advantage of our approach, besides the very high accuracy, is the large range 
of applicability. The convergence of our algorithm does not rely on a good anticipation of 
the bounce and is hence well suited to treat multi dimensional phase transition problems. In 
a forthcoming publication the method we presented will be used to analyze the electroweak 



phase transition in the nMSSM 
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